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We investigate the problem of estimation of parameters of a gravitational 

^^ ■ wave signal from a coalescing binary, at the output of a single interferomet- 

er ' 

ric detector. We present, a computationally viable statistical model of the 

° 

distribution, of the maximum likelihood estimates (MLE), of the parameters. 

This model reproduces the essential features of the Monte Carlo simulations, 
V~j ■ thereby explaining the large root mean square errors in the estimates, ob- 

C3 ■ tained in numerical experiments. We also suggest a more pertinent criterion 

to assess the performance of the MLE, and find that the MLE performs rea- 
sonably. We have considered a Newtonian signal embedded in Gaussian noise, 
with a power spectrum typical of the LIGO/VIRGO type of detectors. The 
model we have used, however, is quite general, and robust, and will be relevant 

to many other parameter estimation problems. 
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Coalescing binaries are relatively clean sources of gravitational waves which can be mod- 
elled with comparatively few parameters. Numerical experiments using simulated signal and 
noise show, that, the actual errors are more than a factor of 3 larger than their Cramer- 
Rao bounds, |1[, at astrophysically relevant signal-to- noise ratios (SNRs). The comparision 
with other less stringent lower bounds [|2[ has also not resolved the discrepancy. In this 
paper we present a statistical model which reproduces the essential features of Monte Carlo 
simulations, even at a 'low' SNR of 7.5. 

The output of the detector will comprise of data trains, each of duration T seconds. We 
assume the number of independent uniformly spaced time samples in each data train to be 
N which is to be chosen such that the Nyquist criterion is met. The set of all detector 
outputs constitutes an N-dimensional vector space V. The statistics of the noise is given by 
the joint probability of the N components of the noise vector n, assumed to be Gaussian, 
and is given by, 



exp 


N/2 

— ~ J2 n J n 3 

. j=~N/2 


l^n 




(2 7 r) iV det 


(-jk 




1/2 



P(n) = / J —'' 1ll/2 J , (1) 



where the components are expressed in the Fourier domain (indicated by a tilde) and Cjj = 
S n (J/T), where S n (f) is the power spectrum of the noise. We assume, S n (f) = (//200)~ + 
2 (l + (//200) 2 ), consistent with the initial LIGO. 

The noise probability distribution introduces a natural metric on the vector space V via 
the scalar product: 

N/2 

(x,y)= £ W 3 IC n . (2) 

j=-N/2 

The set of signal vectors, s(/u) = s(t; /x), where /x = {/i , Hi, • • • , /-tp-i}, is a p-dimensional 
parameter vector, will describe a p-dimensional manifold M. embedded in V. (See |],|3| for 
an introduction to the use of differential geometry in gravitational wave data analysis). Let 
the output of a detector x, contain a signal s(/2). Then x = s(/2) + n, where n is a noise 
vector drawn from the noise distribution. The distance, D(/jl) between x and a point s(/x) 
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is given by, D{(j,) = (x — s(/x), x — s(/^)) ' . The MLE of the parameters is that point p, on 
the parameter space which minimises D(/j,). This is equivalent to maximising the correlation 
c(/x) = (x, s(/x)) provided we keep (s(/x),s(/x)) constant. 

In the limit of high SNR, s(/i) will lie in a small region around s(/2) on the manifold, 
effectively the tangent space to the manifold at that point. In this case, the difference, 
s(/i) — s(/2) can be satisfactorily approximated as a linear function of n. Further, if the 
parameters form a Cartesian system of coordinates, then, they too will be linear in n and 
the distribution of the parameters can be described by a multivariate Gaussian ||. The 
covariance matrix of this distribution defines a lower bound on the errors in estimation and 
is termed as the Cramer-Rao bound. 

If the global minimum of D(fx) is also a local minimum then, at [i = p,, dD([i)/dfi a = 0, 
which implies, 

( S (M)-s( A ),|>) = -(n,|>)). (3) 

Geometrically speaking, ft, is a local extremum when the tip of the vector x lies on that N—p 
dimensional hyperplane Bp,, which passes through s(/t), and is orthogonal to the tangent 
plane at s(/x). This hyperplane is the intersection of the N — 1 dimensional hyperplanes, 
T? each orthogonal to a coordinate basis vector d/dfi a . Let l a be the normalized coordinate 
basis vectors at ft, and let r a be the minimal distance from s(/2) to the hyperplane T? given 
by r a = (sQu) — s(/2), l a ). A schematic illustration of the above is given in Fig. [l|. 

The probability density for the vector x to lie on B^ will depend only on l a and r a , and 
can be written down as, 



p(r a ) 



[p«(n -rj a ),l a » 

a=0 



p{n) d N n. (4) 



Substituting for p(n) in the equation above and integrating, we get, 
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Since l a are tangent vectors on Ai, the metric on the tangent basis is nothing but 
9ab — (la? lb)- We can always select a basis such that g a b are constants over the manifold. If 
the signal manifold is intrinsically flat we can have a coordinate basis with the same property. 
We will assume that the metric coefficients are constant in the \x coordinate system. 

Had the map between ft, to r been bijective, it would have been possible to write the 
distribution for the estimated parameters, simply as, p(r a (ft))J(p,), where, J (p.) is the 
Jacobian of the transformation from fi a to the variables r a . A given set of values for r = {r a } 
could in general correspond to a discrete set of parameter vectors (x^ . The problem now 
is to fix {r a } and compute P(/i^ m ^|r), which is the probability that a particular p,^ will 
be the global maximum for a given r. We will suggest a simple approximation to determine 
P(/i( m )|r), in the context of our problem. 

The gravitational wave from a CB system can be characterised in the Newtonian ap- 
proximation using four parameters, which are, (i) the amplitude of the signal, A, (ii) the 
time of arrival t A , (hi) the initial phase of the the signal (at t A ) <fio, and (iii) the New- 
tonian chirp time r which determines the time left for actual coalesence after t A . Since 
it is convenient to have dimensionless parameters, we define a new set of parameters 
M = {/ib}b=o,...,3 = {A 2irf A t A , 0o, 271-/UT0}, where f A is the frequency at t = t A . An 
expression for the signal in the Fourier domain can be obtained via the stationary phase 
approximation (See for details) as s(f; fjt) = Ah(f;fi k ), with, 
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where, M is a normalization constant such that ||h|| = (h, h) / =1, and ipk(f) are functons 
of frequency as given below: 

Mf) = f/fA, ife(/) = -i, 

*"-H + i(£P (7) 

Henceforth, the indices i,j, k will take values in the range 1 — 3, and the indices a, b in the 
range — 3. 



For the Newtonian chirp, 

1 = h(/i fc ), and l k 
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The coefficients C% b = (l a ,lb) _ are seen to be independent of coordinates for the chirp 
waveform. This is a consequence of the intrinsic flatness of the chirp waveform. In order 
to eliminate the nondiagonal components of C M , we find an orthogonal matrix S which 
diagonalises C M . This transformation leaves the amplitude parameter unchanged. The new 
parameters v are related by the relation v = S/x. In this coordinate system C v ah = 5 a b- 
The root mean square errors as calculated from the covariance matrix, for the parameters 
Uk(k = 1,2,3) at an SNR of 7.5 are {0.020,0.171,4.049} respectively. 
In the v coordinate system, 
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For a fixed set of values, {rj}, we will have multiple solutions, {%} , {pj} 2 \ • • • , {%} , 
to eqn. (W). The correlation obtained at these points will be A"' = (x, h(z>j)^\. (It 
is to be noted that for a fixed {%}, A will be a Gaussian random variable.) We will set, 
P{{^j} \{ r j}) e Q ua l t° the probability that, A® will be larger than the correlation at all 
the other multiple solutions. The identification of the multiple roots is quite a problem, and 
so we make the following approximation. For one of the solutions {%}'-, corresponding to 
{rj}, the probability P({vj}^\{rj}) is almost unity. We assume that {%}^ will be the one 
which is 'closest' to {pj}- If this is true, then, we only need to compute the probability that 
the correlation at an arbitrary point on the parameter space exceeds the correlation at the 
point which shares the same value of {rj} but is 'closest' to {fij}- We obtain this solution 
by linearizing eqn. (|I0| ), to get, 
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So for an arbitrary V}, we follow the following procedure, 
(i) determine rj{i) k ) using eqn. fllPf), 



(ii) determine i> k using eqn. (|ITD , 

(iii) determine the probablility for A = (x, h(r)fc)) to be greater than A® = (x, h(r)fe)^ 

and, 

(iv) Set P({vj}\{tj}) equal to the calculated probability and attach it as a weight functon 

to pip). 

The SNR is essentially the norm of the signal, which is nothing but A. We assume the 
value of 7.5 for the SNR and a detection threshold of 6.5. This effectively means that we 
neglect all events with A < 6.5. Since the amplitude parameter is not of primary interest to 
us, we shall integrate the distribution over A. The distribution p(0k) is finally given as: 
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The integral in the above formula can be evaluated in terms of the error function. 

We now compare the Monte Carlo results with the distribution given in eqn. flT2"l). In 
Fig. ^ we plot the marginal distribution, p(z>2,z> 3 ), obtained by integrating eqn. ( |T2"D with 
respect to v\. The surface plot reveals the locations of the multiple peaks, corresponding 
to the multiple solutions of eqn. ([]]). The periodic and regular spacing between the peaks 
is a consequence of the intrinsic flatness of the intrinsic flatness of the chirp manifold and 
that the selected coordinates are cartesian. We compare the above with the scatter plot in 
Fig. |3] of the estimated values of the parameters on the z/ 2 — u 3 plane, where each point 
corresponds to a realization of noise. The distribution of the plotted points on the z/ 2 — u 3 
plane is consistent with the marginal distribution. Though we have shown only the first 2 
subsidary peaks, the peaks occur even beyond the box shown, but with reduced amplitudes. 
These peaks are again consistent with the Monte Carlo results. 

Each peak in Fig. |^ has a smooth falloff in the z/ 2 direction with an approximately 
Gaussian profile, but there are sharp cutoffs along the z/3 direction. This is due to the 
restriction of the phase parameter to the interval [— 7T, 7r]. The above effect is also seen in 
the scatter plot as the density of points falls of abruptly along the v% direction but not 
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along the 1/2 direction. Further from the inspection of Figures ^| and [^ it is clear that the 
marginal distribution p(z> 2 ) will show pronounced multimodal behaviour, whereas p(0 3 ) will 
not. Moreover, since the parameters fx are linear combinations of u the multimodal nature 
might not be apparent in the fx parameters, as can be observed in |I[. 

We next, compare the one dimensional marginal distribution p{v\) with the histogram 
obtained via the Monte Carlo method. Though the parameter v\ has the least root mean 
square error of .02 as predicted by the covariance matrix, its distribution has the most pro- 
nounced non-Gaussian behaviour. In plot (a) and (b) of Fig. |4| we display the distributions 
p{v\), obtained from the statistical model and Monte Carlo simulations respectively. Plots 
(c) and (d) in the same figure zoom in on the first two maxima occurring after the central 
maximum. It can be seen that in the case of plot (c) the match is not very good even though 
the location of the peaks match fairly well. The discrepancy here could come from either 
the Monte Carlo method or the statistical model. The chief problem seems to be the use of 
the linearized equation to determine v' . Implementing a more accurate inversion procedure 
close to v would be beneficial. 

Since the distribution of the estimated parameters is multimodal, using the variances as 
an indicator of performance of the MLE, is not justified. A more reasonable indicator would 
be to compute the probability that the estimated parameters lie within a certain region 
around the actual value. As a concrete example, we count the number of times (Nr) the 
errors in the estimated parameters in a numerical experiment are less than four times the 
r.m.s. values as computed via the covariance matrix, out of a total, N tot of points. For the 
Newtonian waveform, at an SNR of 7.5 we find that 90% of the points satisfy the above 
criterion. At higher SNRs of 10 and 15 the number of points within the region increases 
to 95.7% and 99% respectively. We also applied this criterion for the results of simulations 
carried out for the first post-Newtonian waveform earlier in |I[ , the results of which are given 
in Table |. 

We see that the MLE works moderately well even at low SNRs. It is to be remarked 
that the assessment of an estimator depends upon how we use the estimates to calculate 



astrophysically interesting quantities. In computing the above probabilities we assume equal 
importance to all parameters. The above gives only a reasonable but general measure of the 
goodness of the MLE. 

We required about 2 days of compution on a 300 MFlops (peak rating) machine to 
generate the results of this paper. The use of an integration routine specifically suited to 
the integrand, and/or the use of lookup tables for computing the integrand, would further 
speed up the computation substantially. Also, for higher values of the SNR the distribution 
will be better behaved, and hence, easier to integrate. The intrinsic flatness of the manifold 
turns out to be very convenient for our purpose. This property holds true for the first post- 
Newtonian waveform as well, and it will not be difficult to implement our model for this case. 
For higher post-Newtonian corrections and/or for inclusion of parameters such as spins, it 
might be computationally expensive to compute the marginal distributions. However, it is 
to be noted that performing Monte Carlo simulations in such cases would also call for a huge 
amount of computational effort. A further research into the above issues is in progress. 
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TABLES 



TABLE I. Nji/Ntot f° r the first post-Newtonian waveform. 



SNR 


10 


12.5 


15 


17.5 


20 


25 


30 


N R /N tot 


0.797 


0.872 


0.922 


0.956 


0.976 


0.980 


0.988 
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FIGURES 

FIG. 1. Schematic illustration of the geometric picture discussed inthe text. 

FIG. 2. Surface plot of the marginal probability distribution, p(^ 2 ,i/ 3 ). The parameter v<i is 
plotted along the X axis, and V3 along the Y axis. The primary peak has been clipped at a value 
of 0.1 at the top. The maximum value attained was 0.185. 

FIG. 3. Scatter diagram of the estimated parameters on the v<i — v% plane. 

FIG. 4. Comparision of the Marginal probability distribution, p{v\) with the histograms 
obtained from the Monte-Carlo simulations. 
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